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Abstract 

We study computational aspects of the nonparametric maximum likelihood estimator (NPMLE) for 
the distribution function of bivariate interval censored data. The computation of the NPMLE consists of 
two steps: a parameter reduction step and an optimization step. In this paper we focus on the reduction 
step. We introduce two new reduction algorithms: the Tree algorithm and the HeightMap algorithm. 
The Tree algorithm is only mentioned briefly. The HeightMap algorithm is discussed in detail and also 
given in pseudo code. It is a very fast and simple algorithm of time complexity O(n^). This is an order 
faster than the best known algorithm thus far, the 0{n'^) algorithm of Bogaerts and Lesaffre (2003). We 
compare our algorithms with the algorithms of Gentleman and Vandal (2001), Song (2001) and Bogaerts 
and Lesaffre (2003), using simulated data. We show that our algorithms, and especially the HeightMap 
algorithm, are significantly faster. Finally, we point out that the HeightMap algorithm can be easily 
generalized to d-dimensional data with d > 2. Such a multivariate version of the HeightMap algorithm 
has time complexity 0{n'^). 
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1 INTRODUCTION 

We consider the nonparametric maximum likelihood estimator (NPMLE) for the distribution function 

of bivariate interval censored data. Let (X, Y) be the variables of interest and let Fq be their joint 

distribution function. Suppose that there is a censoring mechanism, independent of {X,Y), so that 

{X, Y) cannot be observed directly. Thus, instead of a realization (x, j/), we observe a rectangular region 

R(ZV? that is known to contain (x,y). We call R an observation rectangle. Our data consists of n i.i.d. 

observation rectangles Ri, . . . , R„, and our goal is to compute the NPMLE F„ of Fq. 

Let denote the class of all bivariate distribution functions. Furthermore, for each F G JT, let 

PpiRi) denote the probability that the pair of random variables {X,Y) is in observation rectangle Ri 

when {X, Y) ~ F. Then, omitting the part that does not depend on F, we can write the log likelihood 

^Marloes H. Maathuis is a Ph.D. student in the Department of Statistics, University of Washington, Seattle, WA 98195 
(email: marloes@stat. washington.edu) . 
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as 



(1) 



and an NPMLE F„ £ T is defined by 



ln{F„) = max 1„{F). 



As stated here, this is an infinite dimensional optimization problem. However, the number of parameters 



can be reduced by generalizing the reasoning of TurnbuU (19761 for univariate censored data. By doing 
this (see e.g. 



Betensky and Finkelstein (19991, Wong and Yu (19991, Gentleman and Vandal (2001 1) one 



can easily derive that: 

• The NPMLE can only assign mass to a finite number of disjoint rectangles. We denote these 
rectangles by ^i, . . . , Am and call them maximal intersections, following Wong and Yu (1999 1. 



• The NPMLE is indifferent to the distribution of mass within the maximal intersections. 
The second property implies that the NPMLE is non-unique, in the sense that we cannot determine the 



distribution of mass within the maximal intersections. Gentleman and Vandal (2002 I call this represen- 
tational non-uniqueness. Hence, we can at best hope to determine the amount of mass assigned to each 
maximal intersection. Let aj be the mass assigned to maximal intersection Aj, and let a = (qi, . . . , Qm). 
Then PpiRi) in equation ([ij is simply the sum of the probability masses of the maximal intersections 
that are subsets of Ri: 



We can then express the log likelihood in terms of a: 

n n m 

Ua) = ^log(P.(J?0) = ^log(^a,l{^^.c«j). 



(2) 



Let K.^ {ae R'" : Qfj > 0, j = 1, . . . , m} and ^ = {a G : 
ffi™. Then an NPMLE aelCnAis defined by 



l^a — 1}, where 1 is the all-one vector in 



ln(a) = max l„(a). 
aeicnA 



(3) 



This is an m-dimensional convex constrained optimization problem that does not need to have a unique so- 



lution in a. This forms a second source of non-uniqueness in the NPMLE. Gentleman and Vandal (2002 I 
call this mixture non-uniqueness. 

Asymptotic properties of the NPMLE for univariate interval censored data have been studied by 



Groeneboom and Wellner (19921. In contrast to the consistency problems of the NPMLE for bivariate 
right censored data, the NPMLE for bivariate interval censored data has been shown to be consistent 
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(Van der Vaart and Wellner (2000 1). This implies that the NPMLE can be used in practical applications 



to estimate the distribution function of bivariate interval censored data, for example to analyze data from 



AIDS clinical trials (see e.g. Betensky and Finkelstein (1999 1) 



From the discussion above, it follows that the computation of the NPMLE consists of two steps. 
First, in the reduction step, we need to find the maximal intersections A\, . . . , Am- This reduces the di- 
mensionality of the problem. Then, in the optimization step, we need to solve the optimization problem 
defined in 

In this paper we focus on the reduction step. We distinguish between two types of reduction algorithms 
that reflect a trade-off between computation time and space: 

• type 1: The reduction algorithm computes the maximal intersections Ai, . . . , Am- 

• type 2: The reduction algorithm computes the clique matrix, an m x n matrix C with elements 

For n observation rectangles, the number of maximal intersections is O(n^). Hence, given the observation 
rectangles, one can compute the clique matrix from the maximal intersections and vice versa in Oin?) 
time. 

We need 0{n'^) space to store the maximal intersections, while we need O(n^) space to store the 
clique matrix. Thus, type 1 algorithms require an order of magnitude less space to store the output. 
On the other hand, the choice of reduction algorithm determines the amount of computational overhead 
in the optimization step, where the values of the indicator functions l^^^^^.j are needed repeatedly. 
Namely, using a type 1 algorithm requires repeated computation of these indicator functions, while such 
computations are avoided with a type 2 algorithm. Thus, if we use a type 1 reduction algorithm, the 
computational overhead in the optimization step is increased by a constant factor. 

Finally, it should be noted that the clique matrix provides useful information about mixture unique- 
ness of the NPMLE. For example, properties of the clique matrix are used to derive sufficient conditions 



for mixture uniqueness by Gentleman and Geyer (1994 1 and Gentleman and Vandal (2002 I. We can also 



use the clique matrix to describe the equivalence class of solutions to ([Sjl. Let d be a solution, and con- 
sider (C^a) ^ = P&{Rr), i = 1, . . . , n. Since the log likelihood ([2]) is strictly concave in Pa{Ri), the vector 
C"^d is unique. Thus, the equivalence class of NPMLEs is exactly the set {a G /C n .4 : a — C"^&}, 
since all a's in this set yield the same likelihood value. 

We now give a brief overview of existing reduction algorithms. [Betensky and Finkelstein (1999| 



provide a simple, but not very efficient, type 1 algorithm. Gentleman and Vandal (2001 1 introduce a 
type 2 algorithm of time complexity O(n^). Song (20011 proposes a type 1 algorithm that is of com- 



parable speed. The algorithm with the best time complexity so far is the O(n^) type 1 algorithm of 



Bogaerts and Lesaffre (20041. Finally, Lee (19831 gives an 0(71 log n) algorithm for a somewhat differ- 
ent but related problem, namely that of finding the largest number of rectangles having a non-empty 
intersection. 
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In this paper, we introduce two new reduction algorithms. The algorithm we initially developed, 



the Tree algorithm, is only mentioned briefly. It is based on the algorithm of Lee (19831, and is a fast 
but complex type 2 algorithm. Later, we realized the reduction problem could be solved in a much 
simpler way if one is only interested in finding the maximal intersections. This resulted in the HeightMap 
algorithm, a very fast and simple type 1 algorithm of time complexity 0{n^). We discuss this algorithm 
in detail and also give it in pseudo code. Finally, we compare the performance of our algorithms with 



the algorithms of Gentleman and Vandal (20011, Song (20011 and Bogaerts and Lesaffre (20041, using 
simulated data. We show that our algorithms, and especially the HeightMap algorithm, are significantly 
faster. 



2 HEIGHT MAP ALGORITHM 

Recall that we want to find the maximal intersections Ai, . . . , Am of a set of observation rectangles 
-Ri , . . . , Rn ■ There exist several equivalent definitions for the concept of maximal intersection in the liter- 
ature. Wong and Yu (1999 I use the following: Aj 7^ is a maximal intersection if and only if it is a finite 
intersection of the Ri's such that for each i Aj n Ri — 9 or Aj n Ri = Aj. [Gentleman and Vandal (2002[ ) 
use a graph theoretic perspective: maximal intersections are the real representations of maximal cliques 
in the intersection graph of the observation rectangles. 

We view the maximal intersections in yet another way. We define a height map of the observation 
rectangles. This height map is a function /i : ^ N, where h{x, y) is defined to be the number of 
observation rectangles that contain the point {x,y). The concept of the height map is illustrated in 
Figure [1] It is easily seen that the maximal intersections are exactly the local maxima of the height map. 
This is true whenever there are no ties between the observation rectangles, and this observation forms 
the basis of our algorithm. 



2.1 Canonical rectangles 

We represent each rectangle Ri as {x\^i,X2,i,y\,i,y2,i)- The point {x\^i,y\^i) is the lower left corner and 
{x2,i,y2,i) is the upper right corner of the rectangle. We call (a;i,i, 2:2,1] the 3;-interval, and (j/i,i, j/2,i] 
the y-interval of Ri. Furthermore, we use boolean variables (cf i, cf^i, cj* ^, Cj J to indicate whether an 
endpoint is closed. As default we assume that left endpoints are open and right endpoints are closed, so 
that {cl„cl,,cli,cl^) = (0,1,0,1). 

We now transform the observation rectangles Ri, . . . , Rn into canonical rectangles with the same 
intersection structure. We call a set of n rectangles canonical if all x-coordinates are different and all 
{/-coordinates are different, and if they take on values in the set {1, 2, . . . , 2n}. An example of a set of 
canonical rectangles is given in Figure [T] 

We perform this transformation as follows. We consider the a:-coordinates and y-coordinates sepa- 
rately and replace them by their order statistics. The only complication lies in the fact that there might 



4 



column nr 

1 2 3 4 5 6 7 8 9 10 11 12 






1 


1 

R2 


1 


1 


1 




















12 







2 


2 
R3 


2 


1 




















11 







2 


3 


3 


2 


1 


1 

RS 


1 











10 







2 


3 


3 


2 


1 


2 


1 











9 







2 


2 


2 


1 





1 














8 








1 


1 


1 
R4 








1 














7 








1 


1 


2 


1 


1 


2 


1 


1 


1 





6 














1 


1 


1 


2 


1 


1 


1 





5 























1 









R6 





4 























1 








1 


1 


3 
































1 


1 


2 






































1 



u 
SI 



e 

o 



I \ \ \ \ \ 1 \ \ \ 1 \ 1 

123456789 10 11 12 

x-coordinate 

Figure 1: An example of six observation rectangles and their height map. The grey rectangles are the 
maximal intersections. Note that they correspond exactly to the local maxima of the height map. 



be ties in the data. Hence, we need to define how to break ties. We explain the basic idea using the 
examples given in Figure (2] In (a) we have an open left endpoint xi^i and a closed right endpoint X2,j 
with xi^i — X2,j and i ^ j. Then the x-intervals of Ri and Rj have no overlap. Therefore, we sort the 
endpoints so that the corresponding canonical intervals have no overlap, i.e. we let X2j be smaller. In 
(b) we have a closed left endpoint xi^i and a closed right endpoint X2,j with xi^i = X2,j and i ^ j. Now 
the a::-intervals of Ri and Rj do have overlap. Therefore, we sort the endpoints so that the corresponding 
canonical intervals overlap, i.e. we let xi^i be smaller. In this way, we can consider all possible combina- 
tions of endpoints. By listing the results in a table, we found a compact way to code an algorithm for 
comparing endpoints. It is given in pseudo code (Algorithm 1). 

The reason for transforming the observation rectangles into canonical rectangles is twofold. First, 
it forces us in the very beginning to deal with ties and with the fact whether endpoints are open or 
closed. As a consequence, we do not have to account for ties and open or closed endpoints in the actual 
algorithm. Second, it simplifies the reduction algorithm, since the column and row numbers in the height 
map directly correspond to the x- and j/-coordinates of the canonical rectangles. 

2.2 Building the height map 

After transforming the rectangles, we build up the height map. To this end, we use a sweeping technique 
commonly used in the field of computational geometry (|Lee 1983p . By using this technique, we do not 
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Figure 2: Breaking ties during the transformation of observation rectangles into canonical rectangles. 



need to store the entire height map. Instead, we only store one column at a time, in an array hi, ... , h2n. 
To build up the height map, we start with hi, ... , h2n = 0. This is column 1 of the height map. We then 
sweep through the plane, column by column, from left to right. Every time we move to a new column, 
we either enter or leave one observation rectangle. Thus, to compute the values of the height map in 
the next column, we respectively increment or decrement the values in the corresponding cells by 1. For 
example, when we move from the first to the second column of the height map in Figure [T] we enter 
rectangle _Ri. R\ has y-interval (7, 12] which corresponds to rows 8 to 12 in the height map. Hence, we 
increment hg,, . . . ,h\2 by 1. 



2.3 Finding local maxima 

During the process of building up the height map, we can find its local maxima, or equivalently, the 
maximal intersections. We denote the maximal intersections in the same way as the observation rectan- 
gles: Aj = {x\,j,X2,j,yi,j,y2,j). Suppose we apply the sweeping technique to the height map in Figure 
m and suppose we are in column 5. We then are about to leave rectangle R2. The j/- interval of R2 is 
(5, 11], which corresponds to rows 6 to 11 in the height map. Hence, the values of the height map will 
decrease by 1 in rows 6 to 11, and will not change in the remaining rows. Since the values of the height 
map are going to decrease, we may leave areas of local maxima. Therefore, we need to look for local 
maxima in rows 6 to 11 of column 5. We find two local maxima: the cell in row 6, and the cells in rows 
9 and 10. These local maxima in column 5 correspond to local maxima in the height map, say Ai and 
A2 respectively. For Ai, we know that (yi, 1,2/2,1) = (5,6) and for A2 we know that (2/1,2,2/2,2) = (8, 10). 
Furthermore, from the fact that we currently are in column 5, we know that a::2,i = ^2,2 ~ 5. Finally, we 
obtain the values of 2:1,1 and 2:1,2 from the left boundaries of the rectangles that were last entered. For 
the cell in row 6 this is R4 with left boundary 4. Hence, Ai = (4, 5, 5, 6). For the cells in rows 9 and 10, 
we last entered rectangle R3 with left boundary 3. Hence, A2 — (3, 5, 8, 10). From this example we see 
that we need an additional array, ei , . . . , e2n , where ek contains the index of the rectangle that was last 
entered in row k of the height map. 

After finding the first local maxima we can continue the above procedure. However, not every local 
maximum in the array h corresponds to a local maximum in the complete height map. To illustrate this 
problem, suppose that we are in column 6 of the height map in Figure [1] We then are about to leave 
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rectangle Ri with t/- interval (7, 12]. Applying the above procedure, we look for local maxima in rows 8 
to 12 of column 6, and we find a maximum in rows 9 and 10. However, this does not correspond to a 
local maximum in the height map. ft merely is a remainder from the maximal intersection A2 that we 
found earlier. Namely, the local maximum in column 6 is formed by the set {_Ri, i?3} which is a subset of 
the set {Ri,R2, R3} that forms A2. We can prevent the output of such pseudo local maxima as follows. 
After we output a maximal intersection Aj, we set :— for one of the rows in Aj. Then, a local 
maximum in the array h corresponds to a maximal intersection if and only if > for all of its cells. 
In the example in Figure [T] this means that after we output Ai and A2 we need to set et := for one of 
their rows. Ai only consists of row 6, and therefore we set ee := 0. A2 consists of rows 9 and 10, and we 
choose to set eg ~ 0. Then, when we find the local maximum in rows 9 and 10 of column 6, we know it 
does not correspond to a maximal intersection since eg — 0. 

Summarizing, we sweep through the plane from left to right, column by column. At each step in 
the sweeping process we either enter or leave a canonical rectangle. When we enter a rectangle Ri with 
y interval j/2,i], we increment hk by 1 and set :— i for k — yi^i + 1, . . . ,y2,i- When we leave a 
rectangle Ri, we first look for local maxima in hk for k — yi^i + 1, . . . , y2,i. For each local maximum that 
we find in h, we check whether ej, > for all of its cells. If this is the case, we output the corresponding 
maximal intersection and set Ck :~ for one of the cells in the local maximum. Finally, we decrement 
hk by 1 for k — yi^i + 1, . . . ,3/2,1- The complete algorithm is given in pseudo code (Algorithm 2). An 
R-package of the algorithm is available at http://www.stat.washington.edu/marloes. 

2.4 Time and space complexity 

We can easily determine the time and space complexity of the algorithm. In order to transform a set of 
rectangles into canonical rectangles, we need to sort the endpoints of their a;-intervals and y-intervals. 
This takes 0(n log n) time and 0{n) space. At each step in the sweeping process, we need to update 
at most 2n cells of the arrays h and e. Furthermore, we may need to find local maxima in at most 2n 
cells, and we may need to check whether es, > for at most 2n cells. Thus, the time complexity of one 
sweeping step is 0(n). Combining this with the fact that the number of sweeping steps is 0{n) gives a 
total time complexity of O(n^). With respect to the space complexity, we need to store the arrays h and 
e. Hence, the space complexity for computing the maximal intersections is 0{n). However, storing the 
maximal intersections takes O(n^) space. 



3 EVALUATION OF THE ALGORITHMS 



We compared our algorithms with the algorithms of Gentleman and Vandal (20011, Song (20011, and 



Bogaerts and Lesaffre (2004 1, using simulated data. We generated bivariate current status data according 
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Table 1: Mean and standard deviation of the user time in seconds, over 50 simulations per sample size 
from model ([4]). Cells with NA indicate that simulations took over 1,000 seconds to run and were therefore 
omitted. 



to a very simple exponential model: 

X,y,l7,l^~exp(l), (4) 

where X and Y are the variables of interest, U is the observation time for X, V is the observation time 
for Y, and X, Y, U and V are mutually independent. Thus, the observation rectangles were generated 
as follows: 



X, < u„ 


Y, < V, 






= (0,(7„0,yO 


X, < u„ 


Y > V, 








X, > u„ 


Y, < 








X^ > Ur, 






Rr 


= {Ui, oo, Vi,cx} 



We used sample sizes 50, 100, 250, 500, 1,000, 2,500, 5,000 and 10,000. For each sample size, we ran 50 
simulations on a Pentium IV 2.4GHz computer with 512 MB of RAM and we recorded the user times 
of the algorithms. For each algorithm, we omitted sample sizes that took over 1,000 seconds to run. All 
algorithms were implemented in C. 

The results of the simulation are shown in Table [T] We see that the Tree algorithm, and especially 
the HeightMap algorithm are significantly faster than the other algorithms. The HeightMap algorithm 
runs sample sizes of 10,000 in less than two seconds. 

To get an empirical idea of the time complexity of the algorithms, Figure |3] shows a log-log plot of 
the mean user time versus the sample size. We fitted least squares lines through the last 4 points of each 
algorithm. The slopes of these lines can be used as empirical estimates of the time complexity of the 
algorithms. We see that the estimated slope of the HeightMap algorithm is 1.9, which agrees with the 
theoretical time complexity of O(n^) that we derived earlier. Furthermore, we see that the HeightMap 
algorithm is about an order faster than the Tree algorithm, which is about an order faster than the 
algorithm of Bogaerts and Lesaffre. Finally, note that the empirical time complexity of the algorithm of 
Bogaerts and Lesaffre is greater than the theoretical complexity of 0{n'^) that they derived. 
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Comparison of five reduction algorithms 




Figure 3: Log-log plot of the mean user time in seconds versus the sample size, over 50 simulations per 
sample size from model Q. For each algorithm, the estimated slope of its graph is given. These slopes can 
be used as empirical estimates of the time complexity of the algorithms. 

4 MULTIVARIATE HEIGHTMAP ALGORITHM 

The height map algorithm can be easily generalized to higher dimensional data. For example, for 
3-dimensional interval censored data the observation sets Ri take the form of 3-dimensional blocks 
{xi,i, X2,i,yi,i,y2,i, zi^i, Z2,i). In thls situation the height map is a function ft : — > N, where h{x,y,z) 
is the number of observation sets that contain the point {x,y,z). The maximal intersections again cor- 
respond to local maxima of the height map. By first transforming the observation sets into canonical 
sets, this implies that we need to find the local maxima of a 2n x 2n x 2n matrix. We can do this 
by sweeping through the matrix, slice by slice, say along the z-coordinate. We only store one slice of 
the height map at a time, so that h and e are now 2n x 2n matrices. At each step in the sweeping 
process, we either enter or leave an observation set Ri. When we enter an observation set, we update 
the corresponding values of h and e, i.e. we set hk,i := hk.i + 1 and ek.i :— i for all k — xu -f 1, . . . , X2,i 
and I — yi^i -I- 1, . . . ,y2,i- When we leave an observation set, we look for local maxima in the cells of 
the rectangle {xi^i, X2,i,yi,i,y2,i), using the height map algorithm for 2-dimensionaI data. For each local 
maximum that we find, we check whether e^^i > for all of its cells. If this is the case, we output the 
corresponding maximal intersection and set ek^i :— for one of the cells in the local maximum. Finally, 
we decrement hk,i by 1 for k — x\^i 4- 1, ... , X2,i and / = j/i,,; -I- 1, ... , y2,i- 

For d-dimensional data, the time complexity of a sweeping step is 0{n'^~^). Since the number of 
sweeping steps is 0{n), this gives a total time complexity of 0{n'^). With respect to the space com- 
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plexity, we need to store the matrices h and e. Hence, the space complexity to compute the maximal 
intersections is 0{n'^~^). However, storing the maximal intersections takes 0(n'^) space. 



ACKNOWLEDGEMENTS 

This research was partly supported by NSF grant DMS-0203320. The author would like to thank Kris 
Bogaerts, Shuguang Song, and Alain Vandal for providing the code of their algorithms, and a referee for 
suggesting to consider generalizing the HeightMap algorithm to d-dimensional data. Finally, the author 
would like to thank Piet Groeneboom, Steven Schimmel and Jon Wellner for their contributions, support 
and encouragement. 



APPENDIX: PSEUDO CODE 



Algorithm 1: CompareEndpoints(A, 
Input: Two cndpoint descriptors A : 
Output: A boolean value indicating 

Ca:=K,-1) 
CB := (cl, = I) 
rA := {k = 2) 

TB ■■= {I = 2) 

if {xk,i xij) then 

return {xk,i < Xij) 
if {va = Tb and ca = Cb) then 

return {i < j) 
if {va 7^ rB and ca cb) then 

return (r^) 
return (r^ ^ ca) 



B): 

= ixk,z,4,i) and B = ixij,d[j) 
A<B 

{ boolean indicating A is a closed endpoint} 

{ boolean indicating i? is a closed endpoint} 

{ boolean indicating A is a right endpoint} 

{ boolean indicating i? is a right endpoint} 

{ if the endpoints have different coordinates} 

{ . . . then let their coordinates determine their order} 

{ if the endpoints are identical} 

{ . . . then let their index determine their order} 

{ if the endpoints are opposites} 

{ . . . then A < B when ^ is a right endpoint} 

{ otherwise A < B when A is closed left or open right} 
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Algorithm 2: HEiGHTMApALGORlTHM2D(i?i, . . . , i?„): 

Input: A set of n 2-dimcnsional observation rectangles i?i , . . . , i?„ 
Output: The corresponding maximal intersections Ai, . . . , Am 

Transform observation rectangles into canonical rectangles {xi^i, X2.i, yi,i, 2/2, i), using ComparcEndpoints 
Sort xi^i,X2,i, i — 1, . . . , n in ascending order and store their indices i in the list ri, . . . , r2n 
m := { counts number of maximal intersections} 

hi, ... , h2n := { column of height map} 

ei, . . . , e2n ■— Q { index of last entered rectangle; blocks output} 

for j — 1 to 2n do { sweep through height map from column 1 to 2n} 

if {rj is a left boundary) then { we enter a rectangle} 

for k = yi^rj + 1 to y2,rj do { update hk and Cfe for k = yi^rj +!,•••, y2,rj} 

hk ■= hk + 1; ek ■= rj 
else { we leave a rectangle} 

b ■= Vi.xj { bottom coordinate of local maximum; blocks output} 

for k = yi + 1 to y2.rj — 1 do { look for local maxima in rows yi,rj + 1, • • • , 2/2, — 1} 

if {hk+i < hk and 6 > 0) then { there is a local maximum in h} 

if (efi+i, . . . , Cfe > 0) then { the local maximum in /i is a maximal intersection} 

m := m+1; A^ ■= {xi^ekjjjb,k) { output the maximal intersection} 
Cb+i := { set e to zero for row 6 + 1 in A^} 

b:=0 
if {hk+i > hk) then 
b:=k 

k ■= y2,rj { look for local maximum in row 2/2, } 

if {b > 0) then { there is a local maximum in h} 

if (cfo+i, . . . , gfe > 0) then { the local max;imum in /i is a max;imal intersection} 

m:=m+l; Am ■— {xi^ekTjjb,k) { output the maximal intersection} 
Cb+i := { set e to zero for row & + 1 in Am} 

for k = yi^rj + 1 to y2,rj do { update hk for k = yi^^ +!,•••, 2/2,rj} 

hk ■— hk — 1 

Transform the canonical maximal intersections Ai, . . . , Am back to the original coordinates 
return Ai, . . . , Am 
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